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Single cell experiments of simple regulatory networks can markedly differ from cell population 
experiments. Such differences arise from stochastic events in individual cells that are averaged 
out in cell populations. For instance, while individual cells may show sustained oscillations in the 
concentrations of some proteins, such oscillations may appear damped in the population average. 
In this paper we investigate the role of RNA stochastic fluctuations as a leading force to produce 
a sustained excitatory behavior at the single cell level. Opposed to some previous models, we 
build a fully stochastic model of a negative feedback loop that explicitly takes into account the 
RNA stochastic dynamics. We find that messenger RNA random fluctuations can be amplified 
during translation and produce sustained pulses of protein expression. Motivated by the recent 
appreciation of the importance of non-coding regulatory RNAs in post-transcription regulation, we 
also consider the possibility that a regulatory RNA transcript could bind to the messenger RNA and 
repress translation. Our findings show that the regulatory transcript helps reduce gene expression 
variability both at the single cell level and at the cell population level. 

PACS numbers: 87.18.Vf, 87.10.Mn, 87.18. Tt 



I. INTRODUCTION 

The growing interest in biological noise has led to many 
efforts to measure gene expression at the single cell level 
[IH3 , revealing a very distinct dynamics when compared 
to population cell experiments [T]. In two well-studied 
examples, the p53-mdm2 regulatory network and the 
NF-kB signaling pathway, sustained oscillations are ob- 
served in single cells following activation signals [5-7], 
while cell population experiments only show damped os- 
cillations 8-10 . In both cases the core circuit consists of 
a negative feedback loop, one of the most common net- 
work designs, where the active transcription factor pro- 
motes the transcription of its own repressor. 

Stable oscillations are not trivially generated in a single 
negative feedback loop [TT] . A loop composed of only two 
agents does not oscillate for plausible macroscopic equa- 
tions. Sustained oscillations require at least three agents, 
where the third one introduces a time delay that repeat- 
edly causes the system to overshoot or undershoot above 
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and below the steady state [12]. Some models achieve 
sustained oscillations by introducing ad hoc time delays 
to reproduce those that a system incurs when manufac- 
turing the various molecular components [HI 1 13] . The dy- 
namics can also be enriched by considering combinations 
of negative and positive feedback loops [TiHTrJ] , bistable 
switches |17j , or by inheriting oscillatory signals from up- 
stream regulators [lOl [13] . 

In this paper we show that the stochastic fluctuations 
in gene expression in a negative feedback loop can pro- 
duce sustained pulses of protein expression. It has been 
suggested that protein fluctuations are driven by under- 
lying messenger RNA (mRNA) fluctuations [2[ 13] [18]. 
We show that the mRNA stochastic fluctuations can be 
amplified during translation and induce a sustained ex- 
citatory behavior characterized by a series of sustained 
anti-correlated pulses in the expression of the positive 
and negative regulator of the loop. 

Noise induced oscillations have already been found in 
other systems. Oscillations in a circadian clock consist- 
ing of a combination of a positive and a negative feedback 
loop are enhanced by the intrinsic biochemical noise [14] . 
Resonant amplification of the stochastic fluctuations can 
lead to cycling behavior in the Volterra system [TS] and 
in self-regulatory genes [20 . Here we show that a simple 
negative feedback loop consisting of an activator protein 
and its repressor is capable of producing protein pulses 



when the stochastic fluctuations of the mRNA are taken 
into account. This result does not rely on having a large 
number of molecules or on the particular statistical prop- 
erties of the noise, neither depends on upstream pulsating 
signals or couplings to additional loops. 

Recently several experimental studies have shown that 
gene expression occurs in bursts of transcriptional activ- 
ity [HH23]- These bursts are usually ascribed to random 
upstream events, such as chromatin remodeling or ran- 
dom promoter transitions. Here we demonstrate that 
sustained pulses of protein expression can be produced 
merely by the stochastic nature of mRNA kinetics. 

In view of the crucial significance of mRNA fluctu- 
ations, we asked ourselves how gene expression can be 
accurately regulated in the noisy cellular environment. 
Regulation of gene expression is a complex, multi-layered 
process that involves many different players. Since their 
discovery more than a decade ago, regulatory RNAs 
(termed "regRNAs" in this paper) have emerged as key 
regulators in virtually all the cellular processes studied to 
date. regRNAs are non-coding RNAs that regulate gene 
expression by base-pairing to a partially or fully comple- 
mentary mRNA target. MicroRNAs (miRNA) [H] and 
antisense RNA [35] are two examples of regulatory RNAs. 

In the mammalian genomes regulatory RNAs often 
share the same transcriptional regulation as their targets, 
giving rise to a diversity of feed-forward loops [251 I2Z] • 
Such pairs of target and regulatory RNA transcripts are 
found to be co-regulated, co-expressed or inversely ex- 
pressed more frequently than expected by chance, pre- 
sumably due to sharing of common transcription factors 
[351[2?-29 . In particular, some transcription factors have 
been found to bind to overlapping transcript pairs, thus 
potentially coupling the regulation of a coding gene and 
its regulatory RNA |26J. It was suggested that one po- 
tential purpose of such design is to filter transcription 
noise {25) . 

In this work we consider the case in which the posi- 
tive regulator in the loop transcribes both mRNA and 
regRNA. The latter could be either miRNA or an anti- 
sense. Additionally, we assume that the regulatory tran- 
script binds to the mRNA and prevents its translation, 
but without promoting its degradation. We show that 
the presence of regRNA reduces the excitability of the 
system by increasing its capacity to buffer the noise. 



II. MODELS 

We consider three alternative designs of a negative 
feedback loop. The main loop is composed of a transcrip- 
tion factor (the positive regulator, P) and its repressor 
(the negative regulator, N). We assume that in response 
to some external cellular signal, P has been activated and 
is promoting the transcription of N. Two of the designs 
also model the transcript dynamics. In the three models, 
P is degraded at the post-translational level via protein- 
protein interaction with N. Schematic representations 



of the three models, as well as the individual chemical 
reactions, are shown in Fig. [T] The description of the 
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FIG. 1: (Color online) Schematic representation of the re- 
duced, base, and extended models together with the chem- 
ical reactions. The reduced model, making up the simplest 
feedback loop, considers only the production and degradation 
of the proteins (P) and (N), and assumes that mRNA is in 
quasi-equilibrium. The base model includes mRNA dynam- 
ics by expanding the chemical reaction of the production of 
N (red box). The additional chemical reactions describe the 
transcription of sense mRNA (5*), the translation of TV, and 
the degradation of S. Finally, the extended model puts into 
play regulatory RNAs, R (dashed blue box), by adding their 
transcription, degradation, and the formation and degrada- 
tion of a sense-regRNA complex C. 



biochemical parameters and their range of variation are 
summarized in Table II] The values adopted correspond 
to typical mammalian cells. 



A. The base model 

The base model consists of the two proteins regulators 
P and N, and an mRNA transcript S. (The latter will be 
interchangeably refered to as "transcript" , "sense tran- 
script" , or "sense mRNA" .) P promotes the transcription 
of S, which in turn encodes for the negative regulator N. 
The deterministic equations are: 



dP 
~dt 

d A = , 

dN 



= v- a P NP, 



,s, 



dt 



= XS — a N N. 



(la) 
(lb) 
(lc) 



We assume that the cooperative binding of n s molecules 
of P at the promoter site of N are required for efficient 
transcription, and we model the transcription rate with 



TABLE I: Symbols used as parameters of the models, their 
values or range of variation, and their description. 
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Transcription factor (P) induction 
mRNA (S) induction 
anti-sense (R) induction 
protein (A) translation 
sense-regRNA binding 
sense-regRNA unbinding 
A-assisted decay of P 
sense auto-degradation 
regRNA auto-degradation 
protein (A) auto-degradation 
Hill's coeff., sense transcription 
Hill's coeff., regRNA transcription 
Threshold, sense transcription 
Threshold, regRNA transcription 



"dimcnsionlcss parameters 

a Hill's function. The translation rate of N is propor- 
tional to the transcript level, S. Finally, the rate of the 
A^-mediated degradation of P is proportional to the con- 
centration of both proteins. 



rate, and the S-R complex formation and destruction. 
Therefore, we systematically explored how the variation 
in /3 S , /3 R , A, /j on , and fJ, OFF , affect the emergent sys- 
tem bahvior. We maintained the remaining parameters 
constant throughout the simulations; random exploration 
showed that changes in their value did not affect the be- 
havior of the system appreciably. 



C. The reduced model: Neglected RNA kinetics 

Finally, to isolate the effect of RNA fluctuations we 
consider a reduced model where the RNA transcripts are 
assumed to be in quasi-equilibrium. From a molecular 
perspective, this model corresponds to the limit where 
the time scales associated with RNA transcription and 
degradation are much shorter than those associated with 
protein production and degradation. Thus, the tran- 
script levels ad just rapidly to changes in P and N, and 

dS 



tt ~ in Eq. (lb). The equations become: 
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B. The extended model: Incorporation of 
regulatory RNA 

This model introduces a regulatory RNA (regRNA), 
denoted R, that targets S. P promotes the transcription 
of both S and R transcripts. R regulates S by base- 
pairing to it, thus creating a hybrid S-R complex, C. 
We assume that the complex molecule can unbind but 
not degrade, i.e., it returns to the system both comple- 
mentary RNA transcripts. 

The new set of equations describing the model are 
given by: 



dP 
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The extended model contains 14 parameters described in 
Table |l] In this work we focus on analysing how the tem- 
poral behavior of the protein levels is influenced by the 
mRNA and regRNA transcription rates, the translation 



We modeled the stochastic behavior at the single-cell 
level using the Gillespie algorithm [30]. To quantify the 
importance of the stochasticity we solved the determin- 
istic equations and compared their trajectories with the 
stochastic dynamics for each set of chemical parameters 
in all three models. The deterministic equations always 
reached a steady state, which is the same in all three 
models for a given set of parameters. We chose this 
value as initial condition of the stochastic simulations 
in order to minimize the initial transient period. Each 
single stochastic simulation represents a possible cell re- 
alization, the equivalent of a single-cell experiment in 
this description. To simulate the behavior of a random 
distribution of cells, we computed 300 realizations with 
different random seeds for each set of chemical param- 
eters. The typical length of each simulation was 2000 
min. 

For each molecular species X, where X stands for 
P, N, S, R, and, C, we computed the average molecular 
count X = (X(t)) t (averaged across all simulations and 
over time), and the coefficient of variation, V, defined as 
the ratio between the standard deviation and the average 
level 



V(X) 



{[x{t)-xY)\ /2 
x 



(4) 



V significantly increases with the amplitude of pulses 
of X, therefore it provides a general signature of the 



stochasticity and strength of the pulsations for the three 
models. 

To identify the pulsating dynamics in more detail, we 
computed the normalized autocorrelation function |31j . 
given by 



C(r) = 



([X(t)-X][X(t + r)-X]) t 



<[*(*) 



x}\ 



(5) 



The autocorrelation function for a each model and pa- 
rameters was averaged over 300 realizations of the initial 
conditions. 

Additionally, we compared the results of the autocor- 
relation with the power spectrum, defined as S(k) = 
(F(k)F(—k)), where F(k) is the Fourier Transform of 
X(t) and (...) denotes average over realizations. 



IV. RESULTS 

Fig. [2] provides a snapshot of the stochastic behavior 
of the three models. For each model, an example of a 
typical simulation is shown on the left panels, while the 
corresponding correlation functions and power spectra 
(averaged over 300 runs) are depicted on the right pan- 
els. Overall, the base model is characterized by a clear 
excitatory behavior in the positive regulator P. This is 
confirmed by the peaked correlation function, which pro- 
vides a characteristic period of the pulses. At the other 
extreme, the reduced model completely lacks pulses, and 
the correlation function rapidly decays to zero. Finally, 
the extended model shows a more complex behavior: the 
time evolution is characterized by a combination of short 
pulses and long-term fluctuations. The correlation func- 
tion has a smooth, oscillating shape, that corresponds to 
the long-term fluctuations, and a characteristic period of 
the pulses cannot be identified. 

The detailed description of the results for each model, 
as well as an analysis of the conditions for the presence 
or lack of excitations, are provided next. 



A. Stochastic fluctuations of RNAs can produce 
pulses of protein expression 

We first considered the base network, Eqs. ([I]), which 
contains the positive regulator P, the rriRNA S, and the 
negative regulator N. Fig. [2{a) shows the outcome of 
a typical stochastic simulation. For comparison, the P 
deterministic value is also shown. The stochastic and de- 
terministic simulations exhibited remarkable differences. 
While the deterministic concentration was locked in a 
steady state, the stochastic integration showed a sus- 
tained excitatory behavior, where P and N were ex- 
pressed in a series of anti-correlated pulses. The presence 
of regular pulses of protein expression in the base model 
is revealed by the autocorrelation function C(t) of P, as 



shown in the right panel of Fig. [2ja). The autocorrela- 
tion presents a pulsating behavior that provides quantita- 
tive information about the dynamics of the system. The 
negative anticorrelation peak indicates the characteristic 
width of the pulses, around 25 min. The first positive 
peak provides the characteristic interval between consec- 
utive pulses, around 100 min, and that is in agreement 
with the maximum in the power spectrum. 

Minimal transcription and translation rates of the neg- 
ative regulator N were required to identify an excita- 
tory behavior with an observable pattern, as illustrated 
in Fig. [3j Below a minimal values of v and A the amount 
of N was insufficient to implement an effective negative 
feedback: the pulses of protein expression became irregu- 
lar, while the average interval between consecutive pulses 
increased. At extremely low transcription and transla- 
tion rates the pulses disappeared altogether. Above these 
extreme conditions we observed that in general a char- 
acteristic inter-pulse period emerges in the base model 
for all tested sets of parameters (Fig. [3]). The autocor- 
relation shows clear peaks that provide a characteristic 
period. The average period depends on the transcription 
and translation rates, and ranges between 30 and 200 
min for realistic rates. 

The origin of the pulses can be understood by observ- 
ing the mRNA time trace. As shown in Fig. [21(a), the 
mRNA is transcribed in a series of micro pulses that 
are subsequently amplified during translation and inher- 
ited by N. The strong correlation between the pulses of 
mRNA and TV is evident in the figure. 

We also observed a strong correlation between the 
mRNA inter-pulse time lags and the typical width of 
the pulses of P. As seen in Fig. p2[a), P only reached 
significant levels after a substantial decay in N prior to 
a new pulse of N. Therefore the typical length of the P 
pulses depends on the decay time of N, which has been 
taken to be around 5 min [32]. The fast N degradation 
also explains the small width of P pulses compared to 
the typical time gap between consecutive pulses. 

Finally, we note that the excitatory behavior of the net- 
work crucially depends on having a low number of mR- 
NAs. In the base model, the average number of mRNA 
copies was below 5 for most of the explored parameter 
space (Figs. El and [5]). This low average value is in 
agreement with measurements of transcripts copy num- 
ber in mammalian cells, where it was found that many 
transcripts are present in less than one copy per cell on 
average [55] . 



1. Comparison with the reduced model: The importance of 
the mRNA stochastic fluctuations 

To further test the impact of the mRNA stochastic dy- 
namics, we considered an alternative circuit design, the 
reduced model Eqs. (pi), where the mRNA was assumed 
to be in quasi-equilibrium and therefore was not explic- 
itly included in the circuit. Fig. ^[b) shows an example 
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FIG. 2: (Color online) Examples of stochastic simulations in the (a) base, (b) reduced, and (c) extended models. Left panels: 
copy number of P (thick green), N (dotted blue) and mRNA (black). The average P copy number (dashed green line) is 
compared to the deterministic trajectory of P (solid red line). For clarity, their average values are indicated on the right. The 
vertical dashed lines in (a) illustrate the correlation between pulses of N and mRNA. Right panels: detail of the corresponding 
autocorrelation functions C'(r) and power spectra S(k) (inset), averaged over 300 runs. 



of stochastic simulation run with the same biochemical 
parameters as the base model in Fig. [2ta). The pro- 
tein dynamics was considerably different in both models. 
While the pulses were clearly present in the base model, 
no trace of pulses could be found in the reduced model. 
Similar conclusions are reached by comparing the auto- 
correlation and power spectrum in both models. 

Opposed to the base model, intrinsic mRNA fluctua- 
tions were not allowed in the reduced model. Hence, the 
protein levels showed very uniform patterns with only 
small deviations around the deterministic steady state. 
The protein levels were also characterized by a very small 
coefficient of variation. The latter is in agreement with 
the experimental evidence that when the contribution 
of external factor affecting cell-to-cell variability is sub- 



stracted, noise in protein expression is dominated by the 
stochastic production and destruction of mRNAs [21 [3] . 
We stress that both networks were simulated with the 
same algorithm, kinetic parameters, and initial condi- 
tions. 

An additional difference between both models was 
given by the average of the stochastic simulations. The 
average, which provides an insight into cell population 
dynamics, converged towards the deterministic trajec- 
tory in the reduced model, but not in the base model. 
The deviation in the base model became especially im- 
portant at low transcription rates. 

Fig. [4] shows the ratio of the stochastic P copy num- 
ber (averaged over time and over 300 simulations), and 
the deterministic steady state. For transcription rates 
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FIG. 3: (Color online) Color map of the characteristic P 
inter-pulse interval for the base model as a function of the 
transcription and translation coefficients (/3s and A respec- 
tively). A sustained excitatory behavior with a characteris- 
tic period appears once minimum translation and transcrip- 
tion rates are attained. The four surrounding plots, whose 
corresponding (/3 S , A) values are indicated with dots on the 
color map, illustrate the behavior of the autocorrelation func- 
tion C(t), averaged over 300 runs, for 4 extreme combina- 
tions of transcription and translation rates. Except for very 
low transcription and translation rates, the autocorrelation 
C(r) shows peaks of correlation that identify the characteris- 
tic inter-pulse interval. 



below 0.1 transcripts per minute, the deterministic de- 
scription clearly underestimated the real circuit behavior. 
The fact that the divergence appeared at low transcrip- 
tion rates, and that it was absent in the reduced model, 
hints at the mRNA stochastic dynamics as the source of 
this divergence. This is another interesting example of 
a biological network presenting a deviation between the 
stochastic average and the deterministic equations. Such 
deviations are commonly expected in non-linear systems 
and/or systems containing a small number of molecules, 
although more generic systems can also present large de- 
viations 1341. 




Ps (min 1 ) 



FIG. 4: (Color online) Ratio between the P stochastic level 
averaged over time and over 300 simulations and the P deter- 
ministic steady state in the base model, as a function of the 
transcription and translation rates. Contour lines show equal 
stochastic levels of P (black), N (dashed white) and mRNA 
(dotted red). The black cross indicates the location of the 
example shown in Fig. [5Fa) . 



Noise is minimized at high transcription and low 
translation rates 



To achieve a given protein concentration an organism 
can adopt diverse strategies characterized by different 
transcription and translation rates. One strategy con- 
sists of producing a few mRNA transcripts and trans- 
lating them efficiently. Alternatively, the organism can 
transcribe a larger number of transcripts and translate 
each one of them inefficiently. The former strategy is en- 
ergetically favorable since a lower number of transcripts 
has to be produced; however it was suggested that it may 
lead to a noisier pattern of protein expression 35 . 

In order to determine the influence of these two strate- 
gies on noise, we ran simulations for different transcrip- 
tion and translation rates in the base model. For each 
set of parameters we computed the coefficient of varia- 
tion, Eq. Q, as shown in Fig. pjl As a visual help we 
also plotted curves corresponding to equal average levels 
of P, N and S. Examples of stochastic simulations at 
some extreme values are shown in the figure, as well. 

Two axes of variation can be roughly identified. Along 
the main diagonal (from the bottom left to the upper 
right corner) the average level of N grows as the tran- 
scription and translation rates increase. The opposite 
behavior is seen for P, reflecting the negative feedback 
loop between both proteins. At the same time, the coeffi- 
cient of variation of P remains practically constant along 
this axis. Along the other diagonal the coefficient of vari- 
ation of P changes from a minimum at high transcription 
and low translation rates and reaches a maximum at low 
transcription and high translation rates. This behavior 
is seen for the coefficient of variation of N as well (data 
not shown). We also observed that higher levels of noise 
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FIG. 5: (Color online) Color map of the coefficient of varia- 
tion of P in the base model, as a function of the transcription 
and translation coefficients (/3 S and A respectively). Contour 
lines show equal average stochastic levels of P (black), N 
(dashed white) and mRNA (dotted red). P exhibits maxi- 
mum coefficient of variation at low /3 S and high A. The black 
cross indicates the position of the example shown in Fig.[2Ta). 
The four plots illustrate the P(t) behavior for four different 
values of the coefficient of variation, showing that the pul- 
sating behavior of P is maximum at high translation rates 
A. 



are obtained when the average number of mRNA tran- 
scripts is low, and vice versa (see the red equi-(S') lines 
in Fig. [5]), demonstrating the importance of the mRNA 
fluctuations to the excitatory behavior. 

Experimental measurements of protein expression re- 
veal that in many cases the variance in protein levels 
is roughly proportional to the mean. This trend is usu- 
ally explained in terms of an increase in the transcription 
noise with the expression level [2J, [3] . Yet, some genes are 
observed to deviate from this scaling law. 

Our results show a general agreement with the previ- 
ous observation; at a fixed transcription rate, the noise 
increases with the expression level. However, our results 
also show that noise fundamentally depends on the in- 
terplay between the transcription and translation rates. 
For instance, along the lines of constant protein levels the 
coefficient of variation reaches a maximum at low tran- 



scription and high translation rates. This suggests that 
the variability in protein expression depends not only on 
the average expression level, but also on the ratio of tran- 
scription versus translation efficiency. Additional factors 
such as non-linear regulatory loops (as in our model sys- 
tem) may also lead to deviations from the scaling rule. 



B. Regulatory RNA filters transcription noise 

We want to address here the question of how the pres- 
ence of a regulatory RNA transcript, which targets the 
mRNA transcript S, modifies the excitatory behavior in- 
herited from the mRNA stochastic fluctuations. 

Fig. |2jc) shows a typical stochastic simulation of the 
extended model, that includes a regulatory RNA, run 
with the same biochemical parameters as the base model, 
Fig. [2ja) . While pulses of protein expression were present 
in both networks, the extended model, Eqs. (J2|, showed 
broader peaks with significantly lower amplitude relative 
to the average expression. A comparison of the auto- 
correlation in both models showed a weaker spiky behav- 
ior in the extended model. Similarly, while the power 
spectrum shared a similar trend with the base model, it 
lacked the clear maximum present in the latter one. 

The deterministic levels of P (red line) is also shown 
for both models in Fig. [2] While none of the systems 
was correctly described by the deterministic solution, the 
departure from the deterministic dynamics was smaller 
in the model with regRNA. This suggests that, in the 
presence of a regulatory RNA molecule, the individual 
cell's dynamics is more robust to transcript fluctuations, 
hence leading to a reduced cell-to-cell variability. 

To verify that the reduced excitability in the extended 
model is due to the action of the regRNA, we tested 
the dynamics of the extended model at different regRNA 
transcription rates and compared it with the base model 
dynamics. If our hypothesis is correct, the two sys- 
tems should show differences in their excitatory dynam- 
ics as the amount of regRNA increases. Simulations 
started at a minimum sense transcription rate of 0.1 tran- 
scripts/min in order to have a sufficient amount of N for 
effective repression of P. At low regRNA copy numbers, 
the coefficient of variation was maximized at low sense 
transcription rate, in agreement with the qualitative be- 
havior of the base model. However, at high regRNA copy 
number (above 15 transcripts), the coefficient of variation 
became independent of the sense transcription rate. 

The origin of this difference becomes clear when one 
observes the time traces of the different RNA transcripts 
in both models. Fig. [6lb) shows the mRNA time trace 
in the base model, where no regRNA is present in the 
system. Clear mRNA micro pulses with a typical length 
of the order of the mRNA half-life are observed. These 
micro pulses are amplified during translation, producing 
the observed pulses of protein expression. An equiva- 
lent plot for the extended model, Fig. |6jc), reveals a 
completely different dynamics. The regRNA serves as 
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FIG. 6: (Color online) Buffering role of the regulatory RNA in the extended model, (a) Color map of the coefficient of variation 
of P in the extended model, as a function of the sense and regRNA transcription coefficients (/3s and /3 R respectivley) , with 
A = 10. Contour lines show equal average levels of P (black), mRNA (dotted red) and regRNA (dashed yellow). The coefficient 
of variation is maximized at low transcription rates. The black cross indicates the location of the example shown in Fig. [5Ja). 
(b) Example of the mRNA variation (black) for a stochastic simulation of the base model with /3 S = 2.4 and A = 5.0. The 
dotted black line and the dashed red lines show, respectively, the average mRNA concentration and the mRNA deterministic 
trajectory, (c) Top: mRNA concentration (black) for a stochastic simulation of the extended model with identical parameters 
and j3u = 2.0. The mRNA average value is shown with the dotted black line, and the mRNA deterministic trajectory is 
indicated with the dashed red line. Bottom: concentrations of sense mRNA (S, thin black), regRNA (R, dotted brown) and 
S-R complex (C, thick cyan). The right panels show schematic representations of the base and extended models. 



a capacitor: it sequesters sense mRNA when there are 
copies available, and releases them back when there are 
none. Since the typical binding and unbinding rates are 
much faster than the RNA half-life, this sequester and 
release process is repeated many times during a typi- 
cal mRNA micro-pulse, effectively erasing memory from 



any previous mRNA stochastic state. For completeness, 
Fig. [6tc) also shows the concentrations of sense, regRNA 
and sense-regRNA complex. Notice also that the av- 
erage mRNA concentration is significantly closer to the 
deterministic trajectory in the extended model, which ev- 
idences the reduction of the stochastic fluctuations in the 
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network. 

Based on these observations we conclude that the re- 
gRNA is able to filter out some of the stochastic fluctua- 
tions and induce a smoother protein expression pattern. 
Such noise dampening capacity was predicted before in a 
circuit in which a shared transcription factor regulates a 
pair of sense and a highly expressed antisense transcript 

m 

The filtering capability of the regRNA depends cru- 
cially on having fast coupling rates between both RNA 
strands. Fig. f?\ a) shows the coefficient of variation of P 
in the extended model as a function of the sense-regRNA 
binding and unbinding rates, and for two different values 
of the regRNA transcription rate. 

The lowest coefficient of variation is reached in the 
panel with higher regRNA transcription rate. The figure 
also shows curves of equal amount of sense S, regRNA R 
and S-R complex. In both panels the minimum coeffi- 
cient of variation was reached as a trade-off between hav- 
ing fast binding and unbinding dynamics, and having the 
highest number of regRNA molecules (magenta line) . We 
conclude that the noise buffering capability needs both a 
sufficient number of regRNA transcripts and a fast kinet- 
ics to be efficient. This is clearly illustrated in Fig. [Wb), 
which shows examples of stochastic simulations for differ- 
ent combinations of binding and unbinding rates, and for 
two extreme values of the regRNA transcription rate j3 R . 
Indeed, there is practically no buffering for low regRNA 
numbers (left panels), and the time evolution of both P 
and TV shows a pulsating behavior even at large bind- 
ing dynamics. On the contrary, for high regRNA num- 
bers (right panels), such a pulsating behavior is observed 
only at very low binding dynamics, to rapidly disappear 
as soon as the binding-unbinding kinetics, and therefore 
the buffering capacity, increases. 



V. DISCUSSION 

We have analyzed the response of three different net- 
work architectures of a feedback loop and we have ob- 
served remarkable differences among them that emerge 
from the RNA dynamical description. 

The fluctuations inherited from the transcription pro- 
cess can be amplified during translation and produce 
pulses of protein expression. The deterministic approach 
fails to describe such dynamics both at the single cell 
level, where no trace of excitatory behavior is observed, 
and at the cell distribution level, where the deterministic 
steady state underestimates the protein levels. 

Deviations between the stochastic average and the 
deterministic equations may have different origins |34j . 
In our case the difference emerges as a result of the 
non-linear terms, which induce correlations that are not 
present in the deterministic framework. These correla- 
tions become particularly important at low mRNA levels, 
where the largest deviation is observed. Indeed, several 
works in the literature point at the possibility of having 



very few transcripts per cell, in some cases, less than 1 
on average [36) . At such low copy numbers, our find- 
ings stress the importance of using stochastic models to 
accurately describe the network dynamics. 

The excitatory behavior was observed only after simu- 
lating the circuit dynamics with the Gillespie algorithm 
[30] , There are alternative methods to model stochastic 
fluctuations. The Langevin approach, for instance |37j . 
adds a small stochastic term to the continuous determin- 
istic equations to account for noise. This approach is not 
suitable in our system because of the strong RNA fluctu- 
ations, which prevents the characterization of a smooth 
continuous background. An alternative approach is based 
on the linear noise approximation of the master equation 
[38] . It assumes that the system contains a large num- 
ber of particles and models noise as a continuous linear 
gaussian perturbation. This linearized description has 
been applied to processes involving two molecular species 
[39] . Opposed to the linear approaches, the Gillespie al- 
gorithm provides a simple yet powerful method to ob- 
tain stochastic and dynamical solutions compatible with 
the full master equation. It does not require external 
assumptions on the noise, the number of molecules and 
species, or the dynamical regime of the system. 

The protein pulses originate in the RNAs stochastic 
fluctuations. This result is supported by the observed 
correlation between the typical protein pulse length or 
frequency and the mRNA fluctuations, as shown in Fig. [2] 
Additional evidence is provided by the fact that when 
the RNA intrinsic randomness is neglected, the stochas- 
tic behavior converges towards the deterministic concen- 
tration. These results are in agreement with experimen- 
tal evidences that mRNA fluctuations are a fundamental 
source of noise in protein expression [2 [3] . 

To produce a given average protein concentration the 
cellular machinery may choose among different transcrip- 
tion and translation rates. Minimized variability among 
cell population is obtained with a combination of high 
transcription and low translation rates, as verified for 
an auto-regulated gene in steady state using a linearized 
stochastic model [18] . Our results generalize this observa- 
tion. We have analyzed a dynamical, stochastic negative 
feedback loop and found, similarly, that noise is mini- 
mized at high transcription and low translation rates, as 
shown in Fig. [5] 

The presence of a non-coding regulatory RNA may 
help buffer the mRNA fluctuations while allowing the cell 
to maintain low rates of transcription. In this work we 
have considered that a regulatory RNA (regRNA) binds 
to the sense transcript and sequesters it from the cellular 
environment, thus preventing its translation. A regRNA 
molecule sequesters a mRNA when there are some coding 
transcripts available in the medium, and releases them 
back when there are none. If this process is repeated suf- 
ficient times during a typical mRNA half-life, then the 
memory of previous states inherited from the transcrip- 
tion process can be erased and, therefore, noise can be 
partially buffered. In this way, regRNA contributes to re- 
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FIG. 7: (Color online) Effect of the binding-unbinding rates on the filtering capability of the regRNA. (a) Color maps of the 
coefficient of variation of P as a function of the binding and unbinding coefficients (/ioN and /i ff respectively) of the sense 
and regRNA in the extended model. Data are shown for two different regRNA transcription rates /3 R . Sense transcription and 
translation rates are maintained constant. The contour lines show equal average stochastic levels of sense molecules (S, black), 
regRNA (R, dotted pink) and complex (C, dashed white). For /3 R = 0.01 min -1 the number of regRNA is almost constant and 
only the average values is shown, (b) Examples of P (thick green) and N (dotted blue) stochastic simulations for two different 
values of regRNA transcription rate /3r, and for 4 different combinations of complex binding (/Kon) and unbinding rate (/zoff)- 
The sketch in the top-right of each panel indicates the location of the examples shown in (a). The value underneath each 
sketch shows the corresponding coefficient of variation. 



duce the temporal variability at the single-cell level and, 
consequently, also the cell-to-cell variability. A requi- 
site for this mechanism to work efficiently is fast binding 
and unbinding (compared to the typical RNA life-time) 
between the regulatory RNA and its target mRNA, as 
shown in Fig. [6] As a result of this buffering, the ex- 
tended model is less excitable. The circuit still shows 
pulses, yet compared to the base model where there is no 



regRNA, the pulses are broader and with smaller ampli- 
tude, and they appear at a lower frequency. 

A prime example of an oscillating negative feedback 
loop is provided by the p53-mdm2 regulatory network. 
This system has been shown to oscillate at the single 
and cell population levels (sustained oscillations versus 
damped oscillations respectively) 0115]. The maintenance 
and shape of the oscillations have been linked to two up- 
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stream signaling kinases as well as to an additional neg- 
ative feedback loop [10] • Interestingly however, a single 
nucleotide polymorphism (SNP309) found in the mdm2 
promoter that results in higher levels of mdm2 mRNA 
and protein [3U] has been shown to disrupt the oscilla- 
tions of p53 and mdm2 protein [41 . This finding sup- 
ports the idea that RNAs could play an important role in 
the oscillatory dynamics in this network, as a high num- 
ber of mdm2 mRNAs would minimize the importance of 
the intrinsic stochastic fluctuations, and thus attenuate 
the pulses of protein expression. 

While this work has considered a negative feedback 
loop, the mechanism that we have described is more gen- 
eral and apply to a wide variety of regulatory networks. 
We have shown that the RNA dynamics is a fundamen- 
tal source of intrinsic noise, suggesting that a realistic 
description of genetic networks requires the stochastic 



modeling of the transcription stages of protein expres- 
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